Benchmarking GW against exact diagonalization for semi-empirical models 
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We calculate groundstate total energies and single-particle excitation energies of seven pi con- 
jugated molecules described with the semi-empirical Pariser-Parr-Pople (PPP) model using self- 
consistent many-body perturbation theory at the GW level and exact diagonalization. For the total 
energies GW captures around 65% of the groundstate correlation energy. The lowest lying exci- 
tations are overscreened by GW leading to an underestimation of electron affinities and ionization 
potentials by approximately 0.15 eV corresponding to 2.5%. One-shot Go Wo calculations starting 
from Hartree-Fock reduce the screening and improve the low-lying excitation energies. The effect of 
the GW self-energy on the molecular excitation energies is shown to be similar to the inclusion of 
final state relaxations in Hartree-Fock theory. We discuss the break down of the GW approximation 
in systems with short range interactions (Hubbard models) where correlation effects dominate over 
screening/relaxation effects. Finally we illustrate the important role of the derivative discontinuity 
of the true exchange-correlation functional by computing the exact Kohn-Sham levels of benzene. 

PACS numbers: 31.15.bu,33.15.Ry,31.15.V- 



I. INTRODUCTION 

For more than two decades the many-body GW 
approximation of Hedin 1 . has been the state of the 
art for calculating band structures of metals, semi- 
conductors, and insulators 2 ^^. With the entry of 
nanoscience the use of the GW method has been ex- 
tended to low-dimensional systems such as molecules, 
carbon nanotubes, graphene and molecule-surface in- 
terfaces£'£&£ ' 10 i 11 ' 12 . In these systems the interplay 
between quantum confinement (in one or more di- 
mensions) and electronic correlation effects leads to 
novel phenomena like the renormalization of molecu- 
lar electronic levels at surfaces by dynamical polariza- 
tion in the substrat o n i 12 ' 13 i 14 . Very recently, the non- 
equilibrium version of the GW approximation has been 
applied to quantum transport and dynamics in molecular 
junction o 15 ' 16 ! 17 ' 18 ! 19 ' 20 ! 21 ' 22 where dynamic correlations 
seems to be particularly important. 

As the range of systems to which the GW approxima- 
tion is being applied continues to expand, critical inves- 
tigations of the performance of GW for other systems 
than the crystalline solids become important. Here we 
report on benchmark GW calculations for 7r-conjugated 
molecules based on the semi-empirical Pariser-Parr-Pople 
(PPP) mode l 23 ' 24 i 25 . By comparing with exact results we 
obtain a direct and unbiased estimate of the quality of 
the GW approximation in molecular systems. 

Previous benchmark model studies of the GW approx- 
imation have all focused on Hubbard models with local 
interaction s 2 ^ 26 ' 27 ! 28 with the conclusion that GW works 
well for small interaction strengths but fails for larger 
interactions strength. The use of GW in systems with 
local interactions is in fact unfortunate because the im- 
portance of electronic screening, which is the main effect 



described by GW, is weak in comparison to correlation 
effects. In contrast to Hubbard models, the PPP descrip- 
tion includes long range interactions and its parameters 
have been fitted to yield realistic excitation energies of 
conjugated molecules. It therefore provides a better and 
more natural starting point for a study addressing the ac- 
curacy of GW for real molecules and nanostructures. We 
mention that in a related work we have performed first- 
principles GW calculations for a series of 33 molecules 
arriving at very similar conclusions regarding the perfor- 
mance of GW as those reported here^ 

Ab-initio GW calculations typically involve a number 
of "technical" approximations such as the plasmon pole 
approximation, the neglect of off-diagonal matrix ele- 
ments in the GW self-energy, or analytic continuations. 
Moreover they are usually performed non-selfconsistently 
and are subject to basis set errors. In the present work 
the GW calculations are carried out fully self-consistently 
without any further approximations apart from the GW 
approximation itself. 

In this work we calculate total energies and excitation 
spectra of the seven conjugated molecules listed in Tab.fl] 
The excitation spectrum of a system can be obtained 
from the spectral function 
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which has peaks at the excitation energies e n = E n (N + 
1) - E (N) and e n = E (N) - E n (N - 1) correspond- 
ing to electronic addition and removal energies, respec- 
tively. Often in the GW literature, excitation energies 
are often referred to as quasi-particle (QP) energies. In 
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the expressions for the excitation energies E n (N) denotes 
the energy of the nth excited ./V-electron state, \^ n (N)), 
with N referring to the neutral state of the system. For 
molecules the first addition and the first removal energy, 
i.e. n = 0, corresponds to the electron affinity and the 
ionization potential. In Hartrcc-Fock theory Koopman's 
theorem^ states that the eigenvalues of the Hartree-Fock 
Hamiltonian equal the addition/removal energies calcu- 
lated without orbital relaxations in the charged states, 
i.e. e* F = (4$ HF (]V)|ff|4$ HF (iV))-£ HF (JV) for a vir- 
tual orbital n. In particular, the highest occupied molec- 
ular orbital (HOMO) and the lowest unoccupied orbital 
(LUMO) represent well defined approximations to the 
ionization potential and electron affinities, respectively^. 
This approximation neglects two important effects. One 
is the relaxation of the single-particle HF orbitals when 
an electron is removed from or added to the molecule. 
The other is the correlation energy which by definition is 
omitted in HF theory. It is instructive to write the exact 
QP energies as the sum of the three contributions 



A rc l ax "f" A c 



(2) 



The relaxation contribution is the correction that follows 
by calculating the QP energy from self-consistently deter- 
mined HF energies of the neutral and the charged states 
N ± 1. The last term A corr is the remaining contribu- 
tion from the correlation energy. For the addition of an 
electron, i.e. an unoccupied orbital, the relaxation and 
correlation contributions are given by 



A rclax = (N + 1) — Eg* (AO 



(3) 



and 



Acorr = [E n (N + I) - (N + I)] - [Eo(N) - Eq (N)]. 

(4) 

In extended systems the potential due to a single de- 
localized electron/hole decreases with the size of the sys- 
tem. Hence, in such systems there will be no or little 
relaxation of the states due to the addition/removal of 
an electron, and the majority of the correction to the 
QP energy will come from the correlation part A corr . 
In molecules, nanostructures, molecules at surfaces, and 
disordered systems with finite localization lengths, this 
is not the case. Here, the introduction of an additional 
electron or hole will lead to a relaxation of the single- 
particle orbitals corresponding to a screening of the ad- 
ditional charge. As a consequence, the relaxation cor- 
rection A rc i ax to the QP energy cannot be neglected in 
such systems. In fact, we find that A rc i ax is larger than 
A corr for all the molecules studied here, and that the GW 
excitation energies correspond roughly to including only 

Arelax in Eq. ©. 

The paper is outlined as follows. In Sec. |TT] the PPP 
model Hamiltonian for conjugated molecules is intro- 
duced. In Sees. IIII Al and lHI Bl we provide an overview of 
the theory and numerical implementation of the GW and 
exact calculations, and in Sec. IIII CI we discuss the use 



of the von Neumann entropy as a measure of correlation. 
The results for total energies and spectral properties of 
the PPP model are presented in Sees. IIV Al and llVBI and 
a comparison is made to short ranged Hubbard models 
in Sec. IIV CI In Scc lIVDl we calculate the exact Kohn- 
Sham levels for the benzene molecule and compare to the 
exact QP levels. The conclusions are given in Sec. [V] 



II. PARISER-PARR-POPLE HAMILTONIAN 

The Pariser-Parr-Poplc model is an effective ir- 
electron description of conjugated molecules that in- 
cludes electron-electron interactions explicitly. The PPP 
Hamiltonian is given by 



o V v(^i ~ Z i)i^j ~ Z j) + ^2 Uihilflil, (5) 



where cj (ci) creates (annihilates) an electron in the p z 
orbital on atom i of the molecule, hi = + hi{ is the 
number operator, hi a — c\ a c ial Zi is the valence (i.e. the 
number of 7r electrons) of atom i, and (ij) denotes nearest 
neighbour hopping. The Ohno parametrization 3 -^ is used 
for the long range interactions 



14.397 



(28.794/(E/i + U )f 



R 2 



(6) 



where Rij is the inter-atomic distance (in A) and Ui is the 
onsite Coulomb interaction (in eV). For large distances 
the Ohno parametrization recovers the \/r behavior of 
the Coulomb interaction while it for small distances rep- 
resents a screened interaction that interpolates to onsite 
Coulomb interaction Ui for = 0. The onsite energy 
Ej, the hopping clement and the onsite Coulomb in- 
teraction Ui are treated as fitting parameters. In the 
present work values for these parameters have been taken 
from the literatur o 32 ' 33 ! 34 ' 35 ! 36 . Since existing parameters 
have been optimized to optical excitation spectra, an ex- 
act agreement with experimental values for the molecular 
gaps is not to be expected. 



III. METHODS 

A. GW approximation 

Hcdin's equations^ provides a formally exact frame- 
work for the determination of the single-particle Green 
function in a self-consistent manner. In the GW approx- 
imation, which follows by neglecting the so called vertex 
corrections, the electronic self-energy £ is given by the 
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product of the Green function G and the screened inter- 
action W, and can be written symbolically as 

S = iGW, (7) 

where the Green function obeys the usual Dyson equation 
G = Go + GoEG. The screened interaction W is given 
by the bare Coulomb interaction V and the polarization 
in the random-phase approximation (RPA) P = —iGG 
through the Dyson-like equation 

W = V + VPW. (8) 

In fully self-consistent GW the set of coupled equations 
for E, G, P, and W are solved iteratively until the 
Green function has converged. Due to the computa- 
tional requirement of a fully self-consistent GW scheme, 
ab-initio GW calculations are usually carried out non- 
selfconsistently. This approach, which is referred to as 
Go Wo, starts from an approximate Go, typically the non- 
interacting Kohn-Sham Green function, from which a sin- 
gle self-energy iteration is carried out to obtain the final 
Green function. 



1. Numerical details 

The GW calculations have been performed following 
the method described in detail in Ref. |37j. Here we give 
a brief overview of the method for completeness. 

The retarded and advanced single-particle Green func- 
tions are given by 

G r ' a (e) = (e ±irj-H -V H - ^(e))- 1 (9) 

where r\ is a small positive infinitesimal, H$ contains 
the first two terms in Eq. ((SJ), and Vh is the Hartree 
potential. We represent the Green functions and all 
other energy-dependent quantities on a uniform grid, 
—E m , —E m + de, ... , E m . The Fast Fourier Transform 
is used to switch between the energy and time rep- 
resentations. Since 77 determines the minimum width 
of features in the Green function's energy dependence, 
the energy grid spacing should obey de -C 77. All re- 
sults presented here have been converged with respect 
to rj,de,E m . Typical converged values are (in eV) 77 = 
0.02, de = 0.005, E m = 50. 

The lesser/greater Green functions are given by 

G<(e) = -f( £ -p,)[G r -G a ] (10) 
G>(e) = (l-/( £ - /J ))[G''-G Q ] (11) 

where f(s — fi) is the Fermi-Dirac function. The chemical 
potential \i is adjusted to yield the desired number of 
electrons in the system. The formulation in terms of 
a fixed chemical potential rather than a fixed particle 
number is reminiscent of the fact that the method has 
been developed for quantum transport. The one-body 
density matrix is given by 



From p the Hartree and exchange potentials follow 

Vn,ij = 5ij2^2v ik pkk (13) 

k 

V x ,ij = -VijPij, (14) 

where we have defined Vu — Ui, see Eq. ([5]). 

The retarded/ advanced and lesser/greater components 
of the quantities needed to construct the GW self-energy 
read^i 

ScfeiW - iGf^m^it) (15) 
W< f> (e) = £WT fc ( e )P</>( e )Wig.( e ) (16) 

kl 
k 

^ 7> W = WGf (-t) (18) 

The GW equations have been expressed in the time or en- 
ergy domain according to where they are simplest. This 
also reflects the practical implementation. 

The retarded components of Sgw an d P are obtained 
using the fundamental relation 

F r (t) = - l 8(t)lF>(t)-F<(t)} (19) 

which is the Kramers-Kronig relation in the time domain 
relating the imaginary and real parts of F r . 

Since the GW self-energy depends on the Green func- 
tion and vice versa, the equations must be iterated until 
self-consistency. To speed up convergence we use the Pu- 
lay mixing scheme^ as described in Ref. l37l 

2. Total energy 

The total energy can be split into kinetic (and ex- 
ternal), Hartree, and exchange-correlation energy E — 
Eq + En + E xc . In terms of the Green function we have 

E + E K = Tr[H oP ] + -TrfVap] (20) 
For the exchange-correlation energy we have 

E * c = hJ Tr P'» G< ( £ ) + S < (e)G a ( £ )]de, (21) 

where E is the exchange-correlation self-energy. In this 
work E is either the bare exchange, E x , yielding the HF 
approximation, or the GW self-energy, Egw- The expres- 
sion (f2"Tj) follows by expressing (V) in terms of the two- 
particle Green function, G2, and then using the defining 
equation for the self-energy in terms of G^r^- 

B. Exact diagonalization 



(12) 



The most direct way to the spectral properties of a 
system is via the Lehmann representation of the Green 
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function in Eq. ([T]). However, since this requires the full 
set of eigenstates and eigenvalues of the Hamiltonian, 
it is of limited practical use and other routes must be 
taken. The following section gives a brief overview of 
the Lanczos method for iterative diagonalization of large 
matrices. 



1. Calculating the ground state - Lanczos algorithm 

In exact diagonalization the given many-body Hamil- 
tonian is diagonalized directly in the Fock space which is 
spanned by many-particle states (Slater determinants). 
Since the dimensionality of the Fock space grows expo- 
nentially with the number of basis orbitals, symmetries 
of the Hamiltonian can help to reduce the dimensionality 
considerably. For the Pariser-Parr-Pople Hamiltonian in 
Eq. ([5]) the number of up and down electrons, Nt and 
Ni, are good quantum numbers since their correspond- 
ing operators commute with the Hamiltonian. This im- 
plies that the exact diagonalization can be carried out 
in each of the (iVf , iVjJ-subblocks of the Fock space inde- 
pendently. The dimensionality of each (TVj, TVjJ-subblock 
is given by the number of ways iVj spin up electrons and 
Ni spin down electrons can be distributed over L basis 
orbitals, 



(22) 



iV T !(L - JV t )! Ni\{L-Ni)\ 

Very often the ground state is located in the half-filled 
subblock, i.e. iVj = iVj, = L/2 where L is the number of 
basis orbitals. For L — 16 the dimensionality of this sub- 
block is d = 165636900, implying that storing a vector in 
double floating point precision requires ~ 1 Gb of mem- 
ory. With such memory requirements a full diagonaliza- 
tion of the Hamiltonian is of course out of reach. If only 
the ground state is needed, iterative methods can be em- 
ployed. The basic idea of iterative methods is to project 
the Hamiltonian onto the Krylov subspace K, generated 
by repeated applications of H on an arbitrary initial state 
|<M, i-e- 

K = span{|0 o ), H\<f> ), 2 |0 O ), ■ • • , H N \<t>o)}. (23) 

In the Krylov subspace the extreme eigenvalues of the 
Hamiltonian converge fast with respect to the size N of 
the subspace, thus reducing the full diagonalization to 
a manageable diagonalization of a N x N matrix, with 

In the Lanczos algorithm^ the Hamiltonian is pro- 
jected onto a specially constructed orthogonalised Krylov 
basis in which the Hamiltonian has a tridiagonal repre- 
sentation. The basis vectors are generated recursively 
as 



\<f) n +i) = H\4> n ) - a n \4> n ) - b 2 n \4> n -i) 
where the coefficient are given by 



(4>n\4>n 



and b„ 



{4>n-l\4>r, 



(24) 



(25) 
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with initial conditions bo = and \4>-i) — 0. At any 
point during the Lanczos iterations only three Lanczos 
vectors needs to be kept in memory, which makes the al- 
gorithm memory efficient. In the basis of the normalized 
vectors (the basis vectors above are not normalized) the 
Hamiltonian has the following tridiagonal representation 



(26) 



which can be readily diagonalized with methods for tridi- 
agonal matrices. In practice the Lanczos iterations are 
continued until the desired eigenvalues have converged to 
a given tolerance. For the ground state energy Eq, typi- 
cal values for TV range from a few to ~ 200 depending on 
the system size. 

The ground state resulting from a diagonalization of 
the tridiagonal Hamiltonian in Eq. (|26|) is provided in 
the Lanczos basis, i.e. |\&o} = Sn c ™l^™)- ^ n or d er to be 
able to calculate the Green function, its representation 
in the original many-body basis is required. Since the 
Lanczos vectors are not stored, the Lanczos iterations 
must be repeated (starting from the same initial vector) 
to obtain the expansion coefficients cti — ^2 n c n (^i\4>n) 
in the original many-body basis {|$,)}f =1 . 

The most time consuming part of the Lanczos algo- 
rithm is the matrix- vector multiplication H\<f> n ). An ef- 
ficient implementation of this part is hence crucial. For 
this purpose it is convenient to use the bit representation 
of an unsigned integer to code the basis states. Denoting 
the integers with bit representations corresponding to the 
spin up and spin down occupations of a given basis state 
with /j and ij , respectively, the integer representation of 
the basis state is I = J-f + 2 L 7j. With the binary rep- 
resentation of the basis states, the multiplication of the 
Hamiltonian can be done efficiently using bitwise opera- 
tions. 



2. Calculating the Green function 

Having obtained the ground state, the Green function 
can now be calculated. From the Lehmann representa- 
tion it follows that it can be written as 

0^(6)^0^(6) + G%{s) (27) 

with the electron and hole Green functions defined by 

1 



e-H + E$ + it} ll 01 



and 



G%(6) = (^\c) 



6 + H - E£ + it] 



(28) 



(29) 
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respectively. In the following we focus on the electron 
Green function which is the matrix representation of the 
resolvent operator (z — H)^ 1 in the basis spanned by 
the \i) = cI^q) vectors. To obtain the i'th diagonal 
element, 

Gt l {e) = {i\{z-H)- 1 \i), (30) 

where z = e + Eq + irj, again the Lanczos algorithm is 
used to put H on a tridiagonal form, but this time the 
Lanczos iterations are started from the normalized initial 
state \<f>o) = \i)/bo where b\ — Hence, in the gen- 

erated Krylov subspace the diagonal element in Eq. (|30"|) 
corresponds to the matrix element bf\(e — H + Eq + 
ir])~ 1 ]n of a tridiagonal matrix, which can be obtained 
as the continued fraction^ 

G e u (e) = . (31) 

s-a ^_ 

e — a\ 

£ — a<i — ■ ■ ■ 

Again the Lanczos iterations are continued until the fre- 
quency dependent Green function element has converged. 

C. Von Neumann entropy 

The following section demonstrates how a quantitative 
measure of the degree of correlations in a system can be 
obtained by considering the von Neumann entropy of the 
reduced single-particle density matrix p. The entropy is 
defined by 

S[p] = -Tr[plogp] = -5>„logp„, (32) 

n 

where in the last equality p has been expressed in its 
diagonal representation, p = J2 n Pn\ n ){n\- 

In the basis of the atomic p z orbitals the matrix ele- 
ments of the reduced density matrix are given by (with 
the spin index suppressed) 

Pij = (* \ct Ci \* ), (33) 

with the diagonal elements equal to the site occupations. 
In the diagonal representation p n thus represents the oc- 
cupation of the eigenstate \n) of the density matrix. 

We note that < S < L\og2, where 2L is the di- 
mension of the single-particle Hilbert space including 
spin. The expression for S m&x follows because the num- 
ber of electrons equal L in all the systems, i.e. half filled 
"band". When \^q) is a single Slater determinant (cor- 
responding to zero correlation) we have S — 0, and when 
l^o) has equal weight on a complete set of orthogonal 
Slater determinants (corresponding to maximal correla- 
tion) we have p n — 1/2 for all n and thus S = Llog2. 
Thus the number < S/ S'max < 1 represents a natural 
measure of the degree of correlation in l^o)- 



IV. RESULTS 
A. Total energies 

We first address the degree of correlation in the ex- 
act ground states by considering the von Neumann en- 
tropies of the corresponding density matrices. The cal- 
culated entropies are listed in Tab. U Except for the 
Hubbard description of benzene (see Sec. ITVC)) which 
clearly presents strong correlations, the entropies of the 
ground states are ~ 10% of their maximum value S'max 
corresponding to weakly correlated systems. The finite 
values of the entropies reveal that none of the ground 
states are single Slater determinants implying that the 
Hartree-Fock ground state energies will be larger than 
the exact ones. 

We here follow the usual convention and define the 
correlation energy as the part of the total energy not 
included in Hartree-Fock, i.e. 

-Ecorr = ^cxact — Eftp . (34) 

Fig. Q] shows the exact correlation energies of the neu- 
tral molecules together with those obtained by evaluat- 
ing the total energy from Eqs. (|20p and (|2~lj) with the 
self-consistently determined Green function and GW self- 
energy. 

For the series of molecules considered here the correla- 
tion energy constitute less than 0.5% of the total energies. 
Furthermore, as expected it decreases (in absolute size) 
with the number of atoms in the molecule. Clearly, the 
GW approximation performs reasonably well for all the 
molecules capturing on average 66% of the correlation 
energy. 



B. Spectral properties 

For isolated systems such as molecules, true quasi- 
particles resembling single-particle excitations are char- 





Formula 


L 


Sj iS'max 


Sgap (eV) 


thiophene 


C4H4S 


5 


0.07 


11.19 


pyridine 


C5H5N 


6 


0.11 


10.61 


benzene 


C6H6 


6 


0.10 


11.39 


benzene (Hubbard) 






0.50 




biphenyl 


C12H10 


12 


0.10 


9.24 


naphthalene 


C10H8 


10 


0.11 


8.65 


anthracene 


C14H10 


14 


0.12 


7.06 


OPV2 


C14H12 


14 


0.10 


8.30 



TABLE I: Chemical formula, number of p z orbitals (L) in- 
cluded in the PPP model and exact ground state entropies 
(5) for the listed molecules. 




acterized by having a weight close to unity (for non- 
degenerate levels) in the spectral function, i.e. 

Z n = Y,K*n +1 \4K)\ 2 ~l- (35) 

i 

This is equivalent to saying that there exists an orbital 
\v) so that the excited state (^n +1 \ can be written as 
the single-particle excitation cj,^^). In Fig. [2] we show 
the single-particle density of states (DOS), 

D(e)=Y / Me) (36) 

i 

for the OPV2 molecule on a logarithmic scale. The height 
of the peaks reflects the value of Z n (modulo degenera- 
cies). The HF, and in particular, the GW approximation 
reproduce the lowest lying excitations quite well while 
higher excitations are poorly described. All the peaks 
in the HF spectrum have Z n — 1 while GW does shift 
some spectral weight from the main peaks to tiny satellite 
structures (at higher energies than shown on the plot). 
However, the GW satellites do not correspond to features 
in the exact spectrum. This shows that excitations with 
Z„<1, i.e. excitations which do not have single-particle 
character, are not well described by GW whose main ef- 
fect is to improve the position of the HF single-particle 
peaks. 

In the following we consider the lowest lying single- 
particle excitations of the molecules as obtained with 
Hartree-Fock, Go Wo and self-consistent GW. In the 
Go Wo calculations the starting Green function Go is 
taken to be the self-consistently determined Hartree-Fock 
Green function. Fig. Ogives an overview of the calculated 
excitation energies relative to the exact ones. Energies 
corresponding to electron removal and electron addition 
are located on the negative and positive half of the x-axis, 
respectively. From this plot clear trends in the calculated 
excitation energies emerge. 



Within HF the occupied (unoccupied) levels are sys- 
tematically overestimated (underestimated), and the de- 
viation from the exact values worsens for the higher ly- 
ing excitations. A closer inspection of the figure reveals 
a few HF energies at ~ ±5.0 eV and ~ ±5.7 eV that 
more or less coincide with the exact energies. These are 
the HOMO and LUMO levels of the small single-ring 
molecules thiophene, pyridine and benzene. The good 
agreement with the exact levels for these systems does 
not arise because HF gives a correct description of the 
many-body states and their energies. This was already 
clear from the analysis above which showed that the exact 
eigenstates are not single Slater determinants and hence 
the excitation energies in Eq. @ have contributions from 
both A re i ax and A corr . The good agreement must there- 
fore be ascribed to cancellations between the relaxation 
and correlation contribution to the exact energies (this is 
discussed further in connection with Fig. 2]). 

Both the Go Wo and the GW give consistently better 
energies than HF - in particular for the higher lying ex- 
citations where the absolute errors are reduced to less 
than ~0.4 eV as compared to ~1 eV for HF. For the low- 
lying excitations GW slightly overestimates (underesti- 
mates) the occupied (unoccupied) levels corresponding 
to an overcorrection of the HF energies. 

In order to address the relative contributions from 
Areiax and A corr to the excitation energies in Eq. [[2]). 
we plot in Fig. 0] the difference between the exact gaps 
and the gaps obtained from the (i) Hartree-Fock eigen- 
values, (ii) Hartree-Fock total energy differences with 
self-consistent relaxations in the N ± 1 Slater deter- 
minants taken into account, and (iii) the distance be- 
tween the highest occupied and lowest unoccupied peaks 
in the GW spectral function. By using the expression 
for the quasi-particle energies in Eq. 0, the exact gap 
E ga , p = £lumo ~ £homo can be expressed as 

-Egap = EluMO — £ HOMO + ^fclax + ^corr (37) 
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FIG. 3: (Color online) Energy of the 3 highest occupied and 
3 lowest unoccupied molecular orbitals relative to the exact 
values. While Hartree-Fock underestimates the occupied and 
overestimates the unoccupied levels, self-consistent GW shows 
the opposite trends but deviates on average less from the exact 
result. 
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FIG. 4: (Color online) The HOMO-LUMO gap relative to the 
exact values. In addition to the HF and GW single-particle 
energies, the relaxed Hartree-Fock total energy differences, 
E$ F {N + 1) + E$ F (N - 1) - 2E$ F {N) are also shown. The 
excellent results of HF for the three smallest molecules is a 
result of error cancellation between relaxation and correlation 
contributions. 



where A^ x and Ag* are the gap equivalents of the cor- 
responding quantities in Eq. © and £homo/lumo are 
the Hartree-Fock HOMO/LUMO eigenvalues. By defini- 
tion A^^ x is difference between the gaps obtained from 
the HF eigenvalues and relaxed HF total energy differ- 
ences. In Fig. @] this is given by the vertical distance 
between the (blue) squares and circles. The correlation 
contribution Af^ can be read off as the difference be- 
tween the exact gap (dashed horizontal line) and the re- 
laxed HF total energy gap (blue squares). Inclusion of 
relaxation effects clearly reduces the HF gaps consider- 
ably implying that A^ x < 0. This reduction is due to 
the screening from the orbital relaxation which reduces 
the Coulomb interaction with the added hole or electron 
and hence also the gap. 

In contrast to the HF (eigenvalue) gaps for which the 
agreement with the exact gap worsens as a function of 
the size of the molecules, the GW gaps follow more con- 
sistently the same trend and underestimates the exact 
gaps with 0.05 — 0.35 eV for all the molecules. The close 
resemblance between GW and the relaxed HF result in- 
dicates that the effect of GW is mainly to account for the 
screening effects included in HF via orbital relaxations, 

A rc l ax - 



C. Long- versus short-range interactions 

To demonstrate the shortcomings of the GW approx- 
imation for strongly correlated systems, we consider a 
Hubbard model description of the benzene molecule. It 
should be noted that this Hubbard description of benzene 



is not intended as a realistic description of the benzene 
molecule, rather it serves to illustrate the limitations of 
the GW approximation. The Hamiltonian is identical to 
the PPP-Hamiltonian in Eq. except that the long 
range Coulomb interactions in the third term have been 
omitted. The values for the hopping elements and the 
onsite Coulomb interaction are t = 2.539 and U — 10.06, 
respectively. With a £//i-ratio of ~4 this obviously repre- 
sent a strongly correlated system. The latter is reflected 
in the ground state entropy in Tab. [T] which is 50% of its 
maximum value. 

From the calculated total energies we find that the cor- 
relation energy (not included in Fig. [J) constitutes 10% 
of the ground state energy which is a considerably higher 
fraction as for the PPP descriptions of the molecules. 
The GW total energy captures 88% of the correlation 
energy compared to 66% on the average for the PPP de- 
scriptions. However, from an absolute point of view, the 
GW approximation misses the exact ground state energy 
by 0.48 eV. This should be compared to 0.16 eV which 
is the difference between the exact and the GW ground 
state energy for the PPP description of benzene. 

The poor performance of both Hartree-Fock and GW 
for the spectral properties of the Hubbard benzene is il- 
lustrated in Fig. [5] which shows the spectral function as 
calculated with the two methods together with the exact 
one. Both Hartree-Fock and GW severely underestimates 
the position of the LUMO level and completely misses the 
details of the spectrum at higher energies. 

This clearly demonstrates that GW is of limited rele- 
vance when considering systems where correlation effects 
(A corr ) dominates over screening, or relaxation, effects 

(A r elax)- 
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FIG. 5: (Color online) Single-particle DOS for the Hubbard 
description of the benzene molecule (only on-site interactions 
from the PPP model are kept). Note the logarithmic axis. 



tial, is equivalent to the restriction of the Kohn-Sham 
potential V xc (r) in the real-space formulation of DFT to 
a local potential. 

Due to the high symmetry of the benzene molecule 
all sites in the PPP-Hamiltonian are equivalent implying 
that V^ KS has the same value for all sites. Except for a 
constant shift, the eigenvalues of the Kohn-Sham Hamil- 
tonian are therefore given by those of the hopping part 
of the Hamiltonian. The HOMO-LUMO gap calculated 
from the Kohn-Sham eigenvalues is E^. = 5.08 eV which 
is a severe underestimation of the true gap of 11.39 eV. In 
line with previous studie a 41 i 42 we thus conclude that the 
main reason for the discrepancy between KS eigenvalues 
obtained with approximate xc-functionals and the exact 
orbital energies is due to the derivative discontinuity of 
the exact xc-functional. 



V. CONCLUSION 



D. Exact Kohn-Sham orbital energies 

Within density functional theory (DFT) the eigenval- 
ues of the single-particle Kohn-Sham Hamiltonian, i?KSj 
are often interpreted as physical energies. In principle 
the validity of this approximation depends on the size 
of the derivative discontinuity of the true and unknown 
exchange-correlation (xc-) functional— In practice the 
use of semi-local xc-functionals represents an additional 
approximation. It is of general interest to investigate 
to what extent the disrepancy between KS energies and 
true QP energies result from the use of approximate func- 
tional and to what extent this is a property of the ex- 
act functional. Below we compare the exact Kohn-Sham 
spectrum to the exact QP spectrum of the PPP benzene 
molecule using the lattice version of DFT. 

The lattice version of DFT follows by extending the 
fundamental concepts of standard DFT, such as the 
Hohenberg-Kohn theorem and the Kohn-Sham equa- 
tions, to model Hamiltonians as e.g. the PPP- 
Hamiltonian4£ In this reformulation of DFT the site oc- 
cupations rii replaces the continuous electron density n(r) 
as the fundamental variable that determines the ground 
state properties. The lattice version of the single-particle 
Kohn-Sham Hamiltonian is given by the sum of the hop- 
ping terms (the kinetic energy) and a site dependent 
Kohn-Sham potential V^ KS which is constructed to yield 
the correct site occupations of the ground state, 

H KS = - ]T tij clc j<x + J2 ^f S "i- (38) 

For the present purpose the explicit form of the site po- 
tential V^ KS is not important. The fact that the lattice 
version of the Kohn-Sham potential is an onsite poten- 



We have presented calculations for the total energy and 
charged single-particle excitations in seven conjugated 
molecules described by the semi-empirical PPP model 
within fully self-consistent GW and exact diagonaliza- 
tion. The results show that the GW approximation gives 
a consistently good description of both total energies and 
electronic excitations with a slight tendency to overesti- 
mate (underestimate) the position of the latter for occu- 
pied (unoccupied) levels. We have found that the effect 
of the GW self-energy is similar to the inclusion of or- 
bital relaxations in the N ± 1 final states in Hartree-Fock 
theory. On the other hand the contribution to the exci- 
tation energies coming from correlations in the ground- 
and excited states is less well described by GW. This ex- 
plains why GW tend to reduce electron addition/removal 
energies relative to the HF eigenvalues. It was shown 
that GW does not perform well for systems with short 
range interactions (Hubbard models) where correlation 
effects are dominating over screening/relaxation effects. 
Finally it was shown that the exact Kohn-Sham eigenval- 
ues significantly underestimate the true HOMO-LUMO 
gap of a benzene molecule showing the importance of the 
derivative discontinuity of the exact exchange-correlation 
fucntional. 
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